A Study of the Thermodynamic Properties of a Hot, Dense 
Hadron Gas Using an Event Generator 

Hadro-Molecular-Dynamical Calculation 

Nobuo Sasaki 

■ ■■ ■ Physics Department, Hiroshima University, 
Q . Higashi-Hiroshima 739-8526, Japan 
O' 

CN , (Received April 30, 2001 ) 

-)— > ' 
^ ' We investigate the equilibration and the equation of state of a hot hadron gas at finite 

^^ ' baryon density using an event generator which is designed to approximately satisfy detailed 

\ balance at finite temperatures and finite baryon densities of the hadronic scale (80 MeV 

0\ . < r < 170 MeV and 0.157 fm"^ < ub < 0.315 fm"^). Molecular-dynamical simulations 

' were performed for a system of hadrons in a box with periodic boundary conditions. Starting 

CO ' from an initial state composed of only nucleons with a uniform phase-space distribution, the 

^ , evolution takes place through interactions such as collisions, productions and absorptions. 

>0 . The system approaches a stationary state of baryons, mesons and their resonances, which is 

Cn ' characterized by the value of the exponent in the energy distribution common to the different 

CO I particles, i.e., the temperature. After the equilibration, thermodynamic quantities, such as 

ON , the energy density, particle density, entropy and pressure are calculated. Above T ~ m^, 

^^ ■ the obtained equation of state exhibits a significant deviation from the naive mixed free gas 

^^ ' model. Large values of the entropy per baryon are also notable. In our system, the excitation 

^p \ of heavy baryon resonances and meson production are enhanced simultaneously, and the 

^ . increase of the temperature becomes moderate, but a Hagedorn-type artificial temperature 

Oj' saturation does not occur. The pressure exhibits a linear dependence on the energy density. 

D , SI. Introduction 

^<\ • From the viewpoint of microscopic dynamics, a hot, dense hadron gas has many 

\^ . unique properties, e.g., the coexistence of hght relativistic particles and heavy non- 

H \ relativistic particles, the production and absorption of particle-anti-particle pairs, 

■ - - ' and a large number of degrees of freedom of resonance states. Such features are quite 

different from those of ordinary classical molecular dynamical systems, and properties 
of hadron gases in on-equilibrium and off-equilibrium are highly nontrivial and a very 
interesting topic for theoretical study. In particular the equation of state (EOS) 
and transport coefficients of hot, dense hadron gases are quite important quantities 
both in high-energy nuclear physics and cosmology. In the ultra-relativistic heavy- 
ion experiments at CERN and BNL, though the primary purpose is the search for 
a quark-gluon plasma (QGP), the physics of a hadron gas dominates the resulting 
system, and knowledge of the EOS and transport coefficients of a hadron gas is highly 
necessary for a better understanding of the experimental results. ^' In cosmology, 
inhomogeneous nucleosynthesis at an early stage of the universe is considered one of 
the possibilities for the origin of matter, and the baryon diffusion constant is a key 
quantity in this scenario. '^' 

In spite of their great importance, the EOS and transport coefficients of hot, 



dense hadron gases are still poorly known, because of the nonperturbative nature 
of the strong interaction. Numerical simulations based on lattice QCD represent 
a powerful tool for non-perturbative QCD, and many calculations of EOS for hot 
quark-gluon plasma phase have been carried out. ^' ' ^' However, quantitative in- 
vestigation of the hadron phase is very difficult, and only few useful results have 
been obtained. In spite of several novel approaches, ^^' '^' inclusion of the chemical 
potential is still difficult for SU{3) lattice QCD. Moreover, progress in the study of 
transport coefficients is very slow, and only a calculation for the pure gluonic matter 
has been completed. ^' 

Phenomenologically, many thermodynamic models have been proposed. Though 
these models can describe some aspects of the properties of the hadronic matter, 
whether they are realistic enough or not is unclear. For example, Hagedorn proposed 
a bootstrap model many years ago, but the problem of the limiting temperature has 
been pointed out. ^' A hot, dilute hadron gas is often regarded naively as an ideal 
gas of massless pions and massive baryons. Since a pion is the lightest mode to 
be predominantly excited, this picture seems to provide a reasonable starting point. 
However, residual interactions may not be negligible. This simple picture can be 
partly improved by taking the "size" of hadrons into account through an excluded 
volume effect, as in the van der Waals equation. "'^''' ^^^ However, it is still not 
clear how appropriate such an approximation is. Although the excluded volume 
imitates the effect of the repulsion between hadrons, microscopic interactions of a 
hot, dense hadron gas are much more complicated, and their roles are not trivial. 
Thus, we need to investigate the thermodynamic properties of a hadron gas by using 
a microscopic model that includes realistic interactions among hadrons. In this 
work, we have adopted a relativistic collision event generator, URASiMA (ultra- 
relativistic AA collision simulator based on multiple scattering algorithm) , ^^> ' '^^' 
and performed molecular-dynamical simulations for a system of a hadron gas. 

In this paper, we focus our interest on the "hadronic scale" temperature (80 
MeV < T < 170 MeV) and baryon number density (0.157 fm^'^ < ub < 0.315 
fm^^), which are expected to be realized in high energy nuclear collisions. Ther- 
modynamic properties and transport coefficients of hadronic matter in this region 
should play important roles in phenomenological models. Sets of statistical ensem- 
bles are prepared for the system of fixed energy density and baryon number density. 
Using these ensembles, the equation of state is investigated in detail. The statistical 
ensembles have already been applied to calculations of the diffusion constant of a 
hot, dense hadron gas. ^^' Calculations of viscosity and heat conductivity are also 
now in progress. ^^' 

In a previous work, we made a pilot study of the equation of state of a hot, dense 
hadron gas with the event generator URASiMA. ^^' In that work, an approximate 
saturation of the temperature, which behaves like the Hagedorn limiting tempera- 
ture, appears. Recently, Belkacem et al. performed a similar calculation using a 
different event generator, UrQMD, and they also reported the saturation of the tem- 
perature. ^^" '^'^' Though those works have provided valuable information regarding 
the nature of the hadron gas, some results are misleading, because, in those sim- 
ulations, the detailed balance of the processes is explicitly broken. Neither model 



includes multi-body absorptions, which are reverse processes of multi-particle pro- 
duction. For energetic hadrons in a closed system, we need both processes to ensure 
detailed balance. If detailed balance is broken, the reversibility of the equilibrated 
system is no longer realized. A luck of the reversal process of multi-particle pro- 
duction can cause one-way conversion of kinetic energy into particles, and artificial 
saturation of the temperature can occur. Although it is interesting and important to 
formulate these multi-body absorption processes exactly in our simulation, straight- 
forward implementation of them is very difficult and not practical. In this work, 
avoiding this complicated problem, we employed a practical approach. We improved 
URASiMA to almost recover the detailed balance at temperatures of present interest 
*' by adopting an idea that is discussed in the next section. 

In §2, we describe the method of our simulation. We also introduce URASiMA 
and explain the method to improve it so that it almost recovers detailed balance 
effectively. In §3, we discuss the equations of state of the hadron gas, and compare 
them with those of the free gas model. In order to confirm that the obtained result 
is insensitive to the system size, we also investigate the finite size effect. Section 4 
is devoted to a summary and discussion of the outlook. 

§2. Tool and method 



2.1. Method of simulation 

In this work, we focus our investigation on the 
thermodynamic properties of a hadronic system. 
For this purpose, we consider a system in a cu- 
bic box and impose periodic boundary conditions 
in configuration space; i.e., if a particle leaves the 
box, another one with the same momentum enters 
from the opposite side. We display the boundary 
conditions pictorially in Fig. 1, where the box lo- 
cated in the center is the system we refer to as the 
'real box' and the others are replicas. When we 
survey would-be collision points of the particles 
in the real box, we have to survey the collisions 
not only between particles in the real box but also 
with the particles in replicas. 

The energy density e and the baryon number 
density ub in the box are fixed as input parame- 
ters, and these quantities are conserved through- 
out the simulation. The initial distributions of 
nucleons are given by uniform random distribu- 
tions in phase space. The momenta of the particles are adjusted to satisfy the initial 
condition of energy, X]j=i \/^ + pf = e ■ V , va the center-of mass system of the 
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Fig. 1. This picture depicts the 
idea of periodic boundary con- 
ditions. The box that is located 
in the center is the space of in- 
terest. We call this the 'real 
box'. 



*' On the order of several hundred MeV. 



particles: J2i=i Pi — 0- After setting the initial state in the above manner, time 
evolution as described by the event generator URASiMA takes place. Though the 
initial particles are only nucleons, many mesons are produced through interactions. 

In order to confirm the realization of the stationarity of the system, we monitor 
particle densities and collision frequencies. We also check energy distributions of 
each particle. According to these methods of confirmation, as time increases, the 
system seems to become stationary. Whether the energy distributions approach 
the Boltzmann distribution is a fundamental question. If the system is in thermal 
equilibrium, the slope parameters of the energy distributions for all particles should 
converge to the same value, i.e., the inverse of temperature. In order to confirm this, 
we analyze the temporal evolution of the inverse slopes of various particles. 

Running URASiMA many times with the same input parameters and taking 
the stationary configuration in equilibrium, we can obtain statistical ensembles with 
fixed temperature and fixed baryon number density. By using these ensembles, we 
can calculate thermodynamic quantities, such as the particle density, pressure, and 
so on, as functions of temperature and baryon number density. 

2.2. Event generator URASiMA 

URASiMA was originally designed to simulate ultra-relativistic AA collision 
experiments. (Here we call it URASiMA-A.) '^^' For this investigation, we improved 
URASiMA to study thermo-equilibrium systems in a box (URASiMA-B). In both 
models, collisions are realized when the distance between particles becomes smaller 
than their interaction range, which is defined by the relevant total cross-section. 
We describe the fully relativistic method to search for would-be collision points 
in Appendix A.l. The fundamental processes in the URASiMA- A/B are 2-body 
elastic and quasi-elastic collisions between hadrons, multi-particle productions, and 
strong decays of resonances. Hadronic experimental data are used as an input to the 
model. '^^' The processes which URASiMA-B includes are as follows: 

NN ^ NR, (1) 

NN ^ Z\i232/^1232, (2) 

R ^ Ntt, (3) 

R ^ R'tt, (4) 

R ^ Nr, (5) 

r ^^ TTTT, (6) 

NN -^ NN + n secondaries (n > 3). (7) 

Here R and r denote baryon and meson resonances, respectively. The list of particles 
in URASiMA-B is given in Table I. For the calculation of cross sections for quasi- 2- 
body processes, we follow the work of Teis et al. '^^' ' ^^' Cross sections of resonance 
absorptions such as NR -^ NN are calculated using the reciprocity of the S matrix. 
(See Appendix A.2.)^^)"^^) A detailed explanation of multi-particle production is 
given in Appendix A. 3. 



Table I. Baryons, mesons and their resonances included in URASiMA. 



N 


A^938 


iVl440 


Afl520 


A^1535 


A^1650 


iVl675 


A''l680 


A''l720 


A 


-41232 


^1600 


^1620 


^1700 


^1905 


"41910 


"4i95() 




meson 


TX 


0"8()() 


P770 













In the case of high-energy hadron/nuclear colhsions, absorption processes are 
not so important. URASiMA-A includes direct multi-particle production {n > 1), 
and l-vr production/absorption via Zii232- In spite of the limited number of processes 
and particle species, URASiMA-A reproduces the global features of the experimental 
data quite well. ^'^^ 

However, in studies of thermo-equilibrium systems, whether detailed balance is 
maintained or not is an important property of the generator. Though the contribu- 
tions of the multi-particle productions dominate the system at early stages of the 
non-equilibrated system, the reverse process plays an important role in the later, 
equilibration stage. The absence of reverse processes leads to one-way conversion of 
the energy to particles and an artificial temperature saturation in the equilibrated 
system. ^^'' ^^^ However, the exact treatment of multi-particle absorption processes 
is very difficult. In order to treat them effectively with URASiMA-B, direct l-vr 
and 2-7r productions are completely replaced by successive processes of quasi-elastic 
collisions and decays of resonances, and l-vr and 2-7r absorptions through resonances 
are naturally included in the model. For this purpose, URASiMA-B contains A 
and N* particles whose masses are up to 2 GeV. We successfully reproduced I-tt, 
2-7r and inelastic cross sections of NN collisions up to ^/s = 3 GeV without direct 
productions. This is the main difference between URASiMA-A and URASiMA-B. 
For \/s > 3 GeV, in order to obtain an appropriate total cross section, we need to 
take direct production processes into account. In our simulation, only at this point 
is detailed balance broken. Nevertheless, if the temperature is much smaller than 3 
GeV, the influence of the breaking of detailed balance is negligibly small. For exam- 
ple, if the temperature of the system is 100 MeV, the occurrence of such a process is 
suppressed by a factor of exp(— 30), and thus the time scale to detect the violation 
of detailed balance is very much longer than the hadronic scale. 

In order to demonstrate the descriptive ability of URASiMA-B, we compare its 
results with experimental data of the high energy nuclear collisions at BNL/AGS in 
Appendix A. 4. 



§3. Results 



3.1. Parameters 



Focusing our interest on the region of temperature and baryon number density 
which is expected to be produced in high energy nuclear collisions, we used the 
input parameters given in Table II. Here ub = 0.156 fm~^ is taken as the baryon 
number density of normal nuclear matter. Thus, the baryon number density in our 
simulation corresponds to 1 - 2 times larger than that of normal nuclear matter. 
The total isospin of the system is set to zero, i.e., the number of protons is equal to 
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Table II. Input parameters. Here V is the volume of the box, ns is the baryon number density, 
and e is the total energy density. 



V fm"] 


ub [fm ^ 


e GeV/fm^] 




.156 


.250 


.300 .313 .370 


.463 .556 .625 


648 


.741 


.833 .938 


64. 


.234 




.300 


.370 


.463 .556 


648 


.741 


.833 




.313 






.370 


.463 .556 


648 


.741 


.833 




.156 


.250 


.300 


.370 


.463 .556 


648 


.741 


.833 


128. 


.234 




.300 


.370 


.463 .556 


648 


.741 


.833 




.313 






.370 


.463 .556 


648 


.741 


.833 




.157 






.370 


.463 .556 


648 


.741 


.833 


216. 


.231 






.370 


.463 .556 


648 


.741 


.833 




.315 






.370 


.463 .556 


648 


.741 


.833 




.157 


.250 


.300 


.370 


.463 .556 


648 


.741 


.833 


1000. 


.231 




.300 


.370 


.463 .556 


648 


.741 


.833 




.315 






.370 


.463 .556 


648 


.741 


.833 



the number of neutrons at the start. 

In this paper we generated a statistical ensemble through two different methods, 
the phase average (200 event) and the long-time average (1 event). We compared the 
thermodynamic quantities obtained for the two ensembles in the case oi V = 64.0 
fm^, and we confirmed that ergodicity holds to sufficient precision. In this paper, 
thermodynamic quantities are calculated from long-time averages, and the phase 
average is used when we study the time evolution of the system. 

3.2. Chemical equilibration 

Figures 2 (a) and (b) display 
the time evolutions of particle densi- 
ties. These figures show that the sys- 
tem approaches the stationary state 
with time. The saturation of par- 
ticle densities indicates the realiza- 
tion of chemical equilibrium. In or- 
der to confirm the detailed balance 
of each reaction processes, we present 
time evolutions of collision frequen- 
cies for all kinds of reaction processes 
in Fig. 3. As seen them in the later 
stages, the collision frequency of the 
multi-particle production is less than 
10~^ (fm/c)~^, and the time scale 
of this process is much longer than 
that of other processes. Violation of 
reciprocity is important only for such 
long time-scale development. Actu- 
ally, as shown in Fig. 4, no difference 
is observed in the time evolutions of particle densities when the multi-particle pro- 
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Fig. 2. The time evolution of particle densities 
for each particle with V = 64.0 fm^. Two 
results are shown, (a) that for ub ~ 0.156 
\ £ = 0.313 GeV/fm^ and (b) that for 
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Ub 



: 0.156 fm" 
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duction is switched off after t = 150 fm/c. Figure 5 shows that detailed balance 
actually holds on the time scale of several hundred fm/c. We conclude that chemical 
equilibrium in our system is realized. 



n„ = 0.156 fm , e = 0.938 GeV/fm 
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Fig. 3. The time evolution of collision frequencies for various types of collisions: 2-body collisions 
and decays presented in Eqs. (l)-(6), multi-particle productions, and elastic collisions. Here R 
denotes baryonic resonances, and r denotes meson resonances. The calculation was done with 
V = 64.0 fm^, UB = 0.156 fm"^ and e = 0.938 GeV/fm^. 
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Fig. 4. The difference between particle densities for the simulations with and without multi-particle 
production. The calculation without multi-particle production was realized by switching off the 
channel by hand at t = 150 fm/c. The fact that the values in these figures are almost vanishing 
indicates that the contributions of the multi-particle productions are very small. The calculation 
was done with V = 64.0 fm^ ub = 0.156 fm"^ and e = 0.938 GeV/fm^ 
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Fig. 5. The time evolution of the balance of collision frequencies for two different energy densities, 
£ = 0.3f3 GeV/fm"' and e — 0.938 GeV/fm^. In both calculations, the volume and baryon 
number density are V — 64.0 fm"' and ub = 0.156 fm~^. The vertical axes correspond to the 
difference between the frequencies of production processes and absorption processes. The fact 
that the quantities vanish for times greater than t ~ 150 fm/c indicates that the system satisfies 
detailed balance for t > 150 fm/c. 

Here we briefly discuss the time scale of the chemical relaxation of the system. 
Figure 2 shows that the number density n{t) seems to relax exponentially: 



n{t) = {ub — n{oo)) • e ^ + n{oo) for N, 
n{t) = n{oo) • (1 — e~~) for Z\, vr, p. 



(8) 



Based on this fact, we can easily estimate the relaxation time r of particle densities 
that approach chemical equilibrium. We give the results for nucleon, pion and p 
densities in Table III. For the nucleon and pion, the obtained values of r are in the 
range 7-20 fm/c, and they depend on the energy density. The values of r obtained 
for p, however, are significantly larger. 

Table IIL The relaxation time of particle densities for each particle with ub = 0.156 fm~^. 



e [GcV/fm^] 


iV (=Z\) [fm/c] 


■K [fm/c] 


p [fm/c] 


0.313 


7.1 ±0.4 


9.2 ±0.3 


21.2 ±0.7 


0.625 


9.8 ±0.5 


13.3 ±0.3 


44.9 ±1.0 


0.938 


10.6 ±0.6 


18.6 ±0.4 


68.5 ±1.3 



Song and Koch calculated the chemical relaxation time of pions in hot hadron 
gas using the effective chiral Lagrangian. ^^^ They estimated a chemical relaxation 
time for tt at T ~ 150 MeV as 10 fm/c, which is close to the value obtained the 
present results. However, Table HI shows the existence of several different time 
scales, and some of them are very long. Though these values may depend on the 



initial conditions of the simulation, our results indicate the possibility that hadronic 
systems have a long relaxation time in certain cases. 

3.3. Thermal equilibration and temperature 



n„ = 0.156 fm , e = 0.938 GeV/fm' 
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Fig. 6. Energy distributions of Ngss, Zii232, tt and P770 at four different values of time, f = 5 
fm/c, t = 10 fm/c, i = 50 fm/c and t = 200 frn/c. The lines are the fitted results that are 
given by Boltzmann distributions, Cexp(— /3i5). The calculation was done with V = 64.0 fm , 
UB = 0.156 fm~^ and e = 0.938 GeV/fm^ 

Figure 6 displays energy distributions of iVgag, Z\i232, tt and P770 at t = 5, 10, 50 
and 200 fm/c. We plot the results as functions of the kinetic energy K = E — ni, 
so that the horizontal axes for all particle species coincide. The energy distributions 
approach the Boltzmann distribution. 



dN 



dN 



d^p AnEpdE 



Cexp{-pE), 



(9) 



as time increases, where /3 is the slope parameter of the distribution. Moreover, 
the slopes of the energy distributions converge to a common value. These results 
indicate that realization of thermal equilibrium. 

Figure 7 (a) displays the time evolution of the quantities (3~^ that were calculated 
by fitting the energy distributions to a Boltzmann distribution. There, the dotted 
curves correspond to the time evolution of (3^^, i.e., the inverse slope of tt. To 
confirm the establishment of the thermal equilibrium, the difference between the 
inverse slope and that of the pion at time t. 



Ap-\t)=l3-\t)-^-Ht), 



(10) 



was calculated, where j corresponds to Ng^g, A1232 and Pyjq. Figure 7 (b) shows the 
time evolution of these A/3^ . From this figure, it is seen that the Ap^ become zero 
to within the accuracy of the statistics for times later than 150 fm/c. Therefore, we 
conclude that thermal equilibrium is established at about t = 150 fm/c; the values of 
the inverse slope parameters of the energy distribution for all particles become equal 
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for later times. Thus we can regard this value as the temperature of the system. For 
later use, we define the temperature of the system as follows: 



T{nB,e) 






J = ^''938,^i232,vr and ^770, 



(11) 



where the /?• are calculated from energy distributions that are averaged over time 
after t = 200 fm/c, and the Uj denote the errors of the f3~ . 
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Fig. 7. (a) The time evolution of the inverse slopes /3 ^ for Ngsg, /ii232, tt and P770 with V — 64.0 
fm^, UB = 0.156 fm"^, e = 0.938 GeV/fm^. The value of /3"^ was calculated from the fitting 
of energy distributions. Here the dotted curves represent the time evolution of /3^^ for tt. (b) 
The time evolution of Aj3~^ for A^gsg, /ii232 and P770 with V = 64.0 fm^, ub = 0.156 fm^"', 
£ = 0.938 GeV/fm^. Here A(3~^ is defined in Eq. (10). The realization of thermal equilibrium 
can be concluded from these graphs. 



3.4. Thermodynamic quantities 

Figures 8-10 show the relations between the temperature and thermodynamic 
quantities such as energy density, e = ^^^=1^"^ ^'^'^^ Ei, particle density, and pres- 
sure, P = ^J2^=i'^^ '^*^** E • 1^ these figures, all curves correspond to the rela- 
tivistic Fermi-Dirac or Bose-Einstein gas with a finite mass width; i.e.. 



e{T,fi) 
n{T,fi) 
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h 
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where dh is a degeneracy factor and pii{m) is the mass spectral function, which 
is given by the Breit-Wigner distribution for the resonances. In these calculation, 
the baryon chemical potential ^b is adjusted to reproduce the total baryon number 
density, whereas the meson chemical potential pm is fixed to zero. 

In Eq. (12), anti-baryons are ignored, because our event generator does not 
contain anti-baryons. However, their contributions to thermodynamic quantities are 
negligible if the temperature is below 170 MeV, since the ratio of the anti-baryon 
number to the baryon number is suppressed by e t . Even at the smallest value of 
the chemical potential, ^b = 250 MeV (n^ = 0.157 fm"^, T = 170 MeV), the factor 
is about 5% at most. Therefore, quantitative error caused by ignoring anti-baryons 
should be only about several percent. 



Temperature - 



energy density, V = 

0.5 



1 fm 




Figures 8-10 (a) and (b) display 
the equations of state of baryons. In 
these figures, it is difficult to see 
the difference between our results and 
those for the calculation of the free 
gas model. This is the result for the 
when it is under the strict constraint 
of baryon number conservation, which 
fixes the total number densities of the 
A^ and A particles. Thus a close look 
at the fractions of baryon resonances 
is necessary in order to recognize the 
difference between our model and free 
gas model. 

Figure 11 displays the particle ra- 
tios of A''938, Z\i232 and other heavy 
resonances, N* and A*. At lower 
temperature (T = 125 MeV), the dif- 
ference between our results and those 
of the free gas calculations is small. 
However, at higher temperature (T > 
150 MeV), the ratios of light baryons 
(iVgsg and Zii232) become smaller, and 
those of heavy resonances, conversely, 
become larger. Thus we find that the 
influence of interactions clearly ap- 
pears above T ~ m^. Such an en- 
hancement of heavy baryons grows as 

the temperature increases. The maximum value of the enhancement reaches 12-15% 

near T = 170 MeV. 

Figures 8-10 (c) and (d) display the equations of state of mesons. As in the case 

of baryons, the deviation from the free gas model appears at temperatures above 
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Fig. 8. The equation of state of a mixed hadron 
gas at finite temperature (80 MeV < T < 170 
MeV) and finite baryon density (0.157 fm~^ < 
ub < 0.315 fm^ ). Tfie energy densities of (a) 
N, (b) A (tlie sum of all resonances), (c) tt, 
(d) P770 and (e) total are plotted as functions 
of the temperature. The curves correspond to 
the free gas model represented by Eq. (12a). 
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Temperature - particle density, V = 1 fm 
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Temperature - pressure, V = 1 fm 
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Fig. 9. Particle densities of (a) iV, (b) A (the 
sum of all resonances), (c) tt, (d) p77o and 
(e) total are plotted as functions of the tem- 
perature. The curves correspond to the free 
gas model represented by Eq. (12b). 



Fig. 10. Pressures of (a) A'^, (b) A (the sum of 
all resonances), (c) n, (d) prro and (e) to- 
tal are plotted as functions of the tempera- 
ture. The curves correspond to the free gas 
model represented by Eq. (12c). 



T ~ TTT-^. Above this tempera- 
ture, meson production is large, 
and the increase of the temper- 
ature becomes moderate. 

In a previous study, ^^' the 
saturation of the temperature 
appeared, and there we cahed 
it the "Hagedorn-hke hmiting 
temperature." However, it was 
an artificial saturation of the 
temperature, because of the 
luck of the reversal process of 
multi-particle production. In 
the calculation of the improved 
URASiMA, this limiting tem- 
perature does not appear, and 
we believe that this is an im- 
portant result of taking detailed 
balance into account. 



production ratio, V = 10 fm 



T= 125 MeV 



60 
50 
40 
^30 
20 
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n^ = 0.231 fm ■^ 

T= 153 MeV 



• URASiMA 
□ free gas model 



T = 169 MeV 



a 
• 





N.„ A, 



"938 ^1232 




N„. A„„ N A 



938 1232 



938 1332 



Fig. 11. The production ratio of Nqss, "41232 and other 
heavy resonances, N* and A*. Calculations were 
done with ub = 0.231 fm""^ and temperatures T = 
125 MeV, 153 MeV, 169 MeV, respectively. 
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Moreover, our results indicate that the equations of state of baryons and mesons 
are closely related through meson-baryon interactions, such as ttN -^ R and their 
inverse processes. The enhancement of heavy baryon resonances causes the increase 
in the abundances of mesons, and vice versa. Heavy resonances readily produce 27r, 
and thus the enhancement of heavy baryon resonances promotes meson production. 
Therefore, interactions between mesons and baryons are very important in the study 
of the properties of a mixed hadron gas. 

3.5. Entropy 

Figure 12 plots the entropy per baryon versus the temperature. Here the curves 
correspond to free gas calculations. In order to define the entropy, we divide the 
phase space into small cells whose volumes are equal to (7ic)^fm^(GeV/c)^. We 
distinguish each cell by the index /. The density operator of a cell pi = p{x\p^) is 
defined as follows: 



Pi 



Particle exists in the l-ih cell. 
The /-th cell is empty. 



(13) 



The definition of the entropy is given by 

S = -TT{{pi)\n{pi)}, 

all cells in phase space 

= - Yl {pi)Mpi), 

I 
where {pi) is the ensemble average of the density operator of a cell, 

1 



{Pi. 



number of ensemble states 



J2 Pi^ 



(14) 



(15) 



ensemble 



Temperature - S/Ng, V = 10^ fm^ 



40 



20 







• nB = 0.157fm" 
n rig = 0.231 fm ' 
An„ = 0.315fm ' 




80 100 120 140 160 
T [MeV] 

Fig. 12. The entropy per baryon plot- 
ted as a function of the temperature. 
The curves correspond to the free gas 
model. 



and the trace constitutes the average over 
phase space. Figure 12 shows that the 
hadron gas has a larger value of S/Nb in 
the present model than in the free gas model 
for T > m-Tr, because a large part of the en- 
ergy of the system goes into particle (en- 
tropy) production. Moreover, our results 
show a clear dependence on the baryon den- 
sity. Thus, the frequently assumed ansatz 
that S/Nb is insensitive to ub seems un- 
reasonable. Therefore, our results indicate 
that the free gas model provides a poor de- 
scription for these quantities above T = mj^ . 
Hence, we should be more careful in using a 
free gas model in the interpretation of ultra- 
relativistic heavy ion experiments, even if 
some corrections, such as the excluded vol- 
ume effect, are considered. ^^^ " ^'^^ 
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Furthermore, some models seem insufficient in counting the dynamical degrees 
of freedom per baryon. For example, EOS based on the a-u) model ^^-^ predicts 
S/Nb ~ 5 (see Ref. 39)) independently of the temperature, while our result gives 
S/Nb ~ 18.6 ± 0.2 at T = 150 MeV and 29.2 ± 0.1 at T = 170 MeV for the normal 
nuclear matter density. 

3.6. Energy density dependence of the pressure 

From Fig. 13 (a), we find that our results and those of the free gas calculations 
exhibit almost a linear dependence of the pressure on the energy density within the 
range of the present study. Our results can be fitted by the following function with 
parameters h and eo- 

P{e) = h-{e-e^). (16) 



Energy density ■ 



pressure, V 

• Ji, rig 
□ Ji, n^ 
on, n^ 
Ap, rig 
<ap, rig 
Vp, rig 



1 fm 



0.157 [fm 
0.231 [fm 
0.315 [fm 
0.157 [fm 
0.231 [fm , 
0.315 [fm"'; 



0.20 



— 0.15 - 



0.10 - 



Q- 0.05, r 



The values of these parameters 
at three different baryon num- 
ber densities are shown in Table 
IV. 

To clarify the situation fur- 
ther, we now investigate partial 
pressures. As shown in Figs. 8- 
10, pressures and energy densi- 
ties of baryons as functions of 
the temperature do not devi- 
ate from their forms for the free 
gas. Contrastingly, the partial 
pressure of the mesons shows 
a clear deviation from that of 
the free gas model. In Fig. 
13 (b) the partial pressures of 
mesons are plotted as functions 
of the corresponding partial en- 
ergy densities. It is interest- 
ing that pions behave like a free gas with regard to this quantity. On the 
other hand, the slope of the p meson is smaller in the present case than in the 
free gas. This result indicates that the p meson and its interactions play an 
important role in determining the thermodynamic properties of a hadron gas. 



0.00 




0.5 0.7 0.9 
£„ [GeV/fm^] 



0.1 0.2 0.3 

E [GeV/fm'] 



Fig. 13. (a) Pressures are plotted as functions of the en- 
ergy density. Free gas results and a line with slope 
0.21 are also shown, (b) Partial pressures of mesons 
as functions of their partial energy densities. 



3.7. Finite size effect 

Finally, we check the finite size ef- 
fect of our calculation. To see this, 
we prepared four different sizes of the 
boxes, 64.0 fm^, 128.0 fm^, 216.0 fm^ 
and 1000.0 fm"^. The volume depen- 
dence of the total particle density is 
shown in Fig. 14. We were not able 



Table IV. Fitting parameters for Fig. 13. 



UB [fm-^] 


h 


£0 [GcV/fm-^] 


0.157 


0.2171 ± 0.0009 


0.1163 ±0.0026 


0.231 


0.2122 ± 0.0008 


0.1633 ±0.0023 


0.315 


0.2083 ± 0.0007 


0.2188 ± 0.0021 



to find any differences among the four sizes of the box. Other thermodynamic quan- 
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titles yield the same result. As concerns thermodynamic quantities, we can regard 
64.0 fm^ as a sufficiently large spatial size. 



• nB = 0.157 fm"; 
D rig = 0.231 fm"' 
An„ = 0.315fm"' 



E = 0.370 GeV/fm 



E = 0.648 GeV/fm 




2 4 6 8 10 12 
V'" [fm] 



Fig. 14. Volume dependence ol the total particle density for e — 0.370, 0.648 GeV/fm"^ 



We also compared relaxation lengths for various spatial box sizes. In a previ- 
ous paper, we investigated the baryon diffusion constant. ^^' First, the fluctuation- 
dissipation theorem tells us that diffusion constant D is given by the current (veloc- 
ity) correlation,^'^) 



1 r°° 
D = - {v{t)-v{t + t'))dt'. 
o JO 



(17) 



In our definition, the correlation function (• 



(...) 



1 



number of ensemble states 



E 

ensemble 



is given by 
1 



number of baryons 



E 



(18) 



baryons 

The correlation functions are damped exponentially with time (see Fig. 15): 

t' 



(v(t)-v(t + t')) ocexp 



Thus the diffusion constant can be rewritten in the simple form 



D = -{v{t)-v{t))TB, 



(19) 



(20) 



where tb is the relaxation time of the baryon current. Figure 16 shows the volume 
dependence of tb • {vb), where (vb) denotes the average speed of baryons, and thus 
Tb • {vb) has dimensions of length. Though this quantity has a clear dependence 
on the size of the box at low energy density {e = 0.370 GeV/fm^), this box size 
dependence seems to disappear for volumes larger than about 6^ fm^. In this work, 
to be safe, we used a box of volume V = 1000 fm^. Therefore, the result in this 
paper can be considered free of any box-size effect. 
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Fig. 15. Velocity correlation of the baryons as a function of time. The normalization of the data is 
arbitrary. 
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An„ = 0.315fm"' 



8 = 0.370 GeV/fm 



E = 0.648 GeV/fm 




V [fm] 



2 4 6 8 10 12 
V'" [fm] 



Fig. 16. The relaxation time nmltiplied by average velocity as a function of the volume of the 
system. 
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§4. Summary 

In this paper, we studied the equilibration and the equation of state of hot (80 
MeV < T < 170 MeV) hadron gas with finite baryon number density (0.157 fm~^ 
< ub < 0.315 fm~^) by using the event generator URASiMA, which maintains 
detailed balance in the practical sense in the hadronic energy scale. We performed 
molecular dynamical calculations for a system of hadrons in a box with periodic 
boundary conditions. The energy density and baryon number in our simulation 
correspond to those of a hot hadron fluid, which is believed to be produced in high 
energy nuclear collisions. Our results for thermodynamic quantities can be related 
to experimental data through statistical models and hydrodynamic models. ^^' 

Collision frequencies and particle densities exhibit saturation as the time evolu- 
tion proceeds. These results indicate that the system reaches chemical equilibration. 
We also studied the thermal equilibration of the system from the time evolution of 
the inverse slopes of the energy distributions. To confirm that thermal equilibrium is 
established, it was demonstrated that the slope parameters of Ng^s, Z\i232, ^r and pjjo 
become almost identical for times greater than t ~ 150 fm/c. Thus the temperature 
of the system can be defined after this time. 

After the equilibration, thermodynamic quantities were evaluated, and the equa- 
tion of state was investigated. Energy densities, number densities, pressures and 
entropies per baryon were plotted as functions of the temperature. Deviations from 
the free gas model were manifestly observed above T ~ m-,^. Above this temper- 
ature, the excitation of heavy baryon resonances were found to be enhanced, and 
meson production was found to be significant. Those effects suppress the increase 
of the temperature, but the saturation of the temperature, as in the case of Hage- 
dorn's limiting temperature, never occurs. These notable differences from previous 
works ^^^'^^^ are important results of the proper maintenance detailed balance be- 
tween production processes and absorption processes. For a reliable simulation of a 
statistical system, detailed balance is essentially important. In our study, we found 
large values of the entropy per baryon. This depends strongly on the baryon density 
above T ~ m^. The pressures exhibit linear dependences on the energy densities 
within the range of present study. Such behavior was analyzed in detail by looking at 
relations between partial pressures and partial energy densities of mesons. We find 
that the p meson and its interactions play an important role in the thermodynamic 
properties of the hadron gas. 

Because the temperature in this investigation (80-170 MeV) is much lower than 
the mass of a KK pair, the hadron gas is mainly composed of non-strange particles. 
For this reason, strange particles are not considered in our model. However, in 
the early stages of the AA collision, the energy density can be very high, and the 
strangeness degree of freedom should play an important role. Including strange 
particles in our model is the next task, and we are now in process of doing so. 

The investigation reported in this paper was focused on "hadronic scale" energy 
densities and baryon number densities. However, it would be interesting to perform 
calculations at higher energy densities and baryon number densities. Our EOS with- 
out a phase transition can work as a helpful reference to the equation of state with 



a QGP phase transition. 
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Appendix A 

Event Generator URASiMA 



A.l. The search for would-be collision points 

After a collision at the space-time point Xg ^ 
i-th particle is given by the straight line 

fi 

fi /, \ ^i t \_ ^ 

Uli 



{toi,XQi), the trajectory of the 



(Al) 



where pf = {Ei,p^) and mj are the momentum and the mass of the particle after 
the collision, and ^i represents the proper time: 

Ci = -E^it^-to,). (A-2) 



Equation (A-1) holds until the next interaction 
occurs. 

We adopt a hard sphere approximation for bi- 
nary collisions; that is, a collision occurs whenever 
the condition 

Trb*^<a (A-3) 

is satisfied, where |6*| is the impact parameter in 
the CM frame of two particles, and a is the total 
cross section. For the manifestly relativistic for- 
mulation, it is important to express the collision 
condition (A-3) in a covariant way. For this pur- 
pose, we consider the CM frame of two particles, 
1 and 2 (see Fig. 17). When two particles are at 
the points of their closest approach, the impact 
parameter b* = xl{tl) — 3^2(^2) ^^"^ the time t* 
are defined as 



{xuti)-xm)) 

t*l =t*2 = 



■- t* 



0, (i 



1,2) 
(A-4) 



-^n/ff/) J 



-f 



?-4 



4'2^ 



Xn?(^n?' 



where 'j — ■ ^2 

Fig. 17. The condition for which 
the collision occurs in the CM 
frame of two particles. In this 
frame, the momentum vector 
Pi II P2 and the difference be- 
tween positions Xi{tl) — 0:2(^2 ) 
become perpendicular at the 
point of nearest approach. In 
this situation, a collision occurs 
if Eq. (A-3) is satisfied. 
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where t* denotes the time at which a collision may occur. Here, we define the total 
momentum Q^* and the diff'erence of momenta K^* as 

K^* ^ pT - pT = (El - E*„pl - pl). (A-5) 

Using Qf^* and K^^*, Eq. (A-4) can be rewritten as 

{t*i - t*2)Q°* = 0. (A-6) 

The invariant expressions of Eq. (A-6) are easily obtained as follows: 

(xi - X2) -Q = 0, 

(xi -X2)-K = 0. (A-7) 

Replacing xi and X2 in Eq. (A-7) by their forms in Eq. (A-1), we find the equations 



mi 771,2 

Pi- K P2- K 



C2 + Xo-K = 0, (A- 



nil "1-2 

where Xq = Xq^^ — Xq2- The solutions of these equations are easily obtained as 

(|)^5^(:S |)(f2)' <«' 

where J is defined as 

{pi ■ Q){P2 ■ K) - {pi ■ K){p2 ■ Q) 



J 



17111712 
, {Pl ■P2f' - {17111712)'^ 

mim2 



(A-10) 



Using Eqs. (A-1) and (A-9), the invariant expression for the impact parameter is 
given by 

b" = (Xl - X2f 



2 (^0 • Q? [(Aq -k)- q, \ 



This expression enables us to specify the collision point in terms of momenta and 
space-time coordinates of the starting points (the previous collision points). In the 
simulation, all candidates for collision are searched for using Eq. (A- 11), and the 
earliest collision in the rest frame of the box is generated. 
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A. 2. Absorption cross sections 

By using the reciprocity of the S matrix, the absorption cross section can be 
related to the production cross section 



1 



0'N"R^NN' 



9rI + ^nn' pi 



p) 

jaNN'-^N"R, 



(A-12) 



with 6iy]\fi the Kronecker 6, which is equal to 1 when the final state nucleons belong 
to the same state. The relation Eq. (A-12) is limited to particles with definite 
masses. Extension of Eq. (A-12) to take into account the width of resonance mass 
is straightforward, '^^' ' ^'^' and we have 



(^N/'R^NN' 



ruRp'j 



<^NN'^N"R 



5R 1 + ^NN' „. fV^-rriN ^^IJhlm' nnim' ^n 
P^ Jmn+mK 27r ^RPRV^R)Pi 



(A-13) 



where puimj^) is the mass distribution function, which is given by the Breit-Wigner 
distribution. 



A. 3. Multi-particle productions 



Secondary Particles 



<E>= 
a/(l + a 




Fig. 18. The diagram of the multi-particle 
production process for an A^ + A'' in- 
elastic collision. Here (E) denotes the 
average of the energy fraction that the 
primary nucleon carries after the colli- 
sion. This quantity is related to the 
value of a, which is one of the ad- 
justable parameters. The remainder 
of the energy fraction is consumed by 
the production of secondary particles, 
whose momenta are determined ac- 
cording to the distribution character- 
ized by /3. 



In high energy hadron-hadron colli- 
sions, n-particle (n > 3) productions may 
occur. In our model, such multi-particle 
productions are treated as direct processes. 
Such processes are very important in the 
early stages, where energetic particles dom- 
inate the system. In the present simulation, 
initial states include many energetic nu- 
cleons, and multi-particle production takes 
place frequently. In the later stages, how- 
ever, they seldom occur, and their effect is 
negligibly small. 

In URASiMA, a multi-particle produc- 
tion process is realized by use of the multi- 
chain model, ^^-^ " ^^' where the exchange of 
a hadronic chain causes the direct emission 
of multipions. Figure 18 displays a dia- 
gram of A^ -|- A^ inelastic collisions. In this 
figure, nucleons exchange a chain, which 
produces secondary particles. Such multi- 
particle production is specified by param- 
eters as follows. The probability distribu- 
tion of a light-like longitudinal momentum 
fraction xn that is carried by nucleons after 
collisions is given by 



P{xn) 



a ■ X 



a-l 

N ■ 



(A- 14) 
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The average energy that is carried by nucleons is a/(a + 1) times the initial energy. 
The remainder, which is on average equal to l/(a+l), is consumed by the production 
of secondary particles. 

In the CM. frame of produced par- 
ticles, the longitudinal momentum of Table V. Parameters of multi-particle produc- 
the secondary particle behaves accord- 
ing to the distribution 



tion. 



Q 


P 


TO [fm/c] 


1 


1.0 


1.0 



dN ( 2mTCOshy\'3 

with rax and y the transverse mass and the rapidity of the secondary particle. The 
quantity v^ stands for the energy deposited in the blob of the diagram. 

Transverse momentum distributions of primary and secondary particles are given 
by gamma distributions: 

-—^(xpTe~ ^^^- (j = nucleon, vr, etc.) (A-16) 

dpT 

Here pT denotes the transverse momentum, and the slope parameter Bj is another 
parameter of the model. 

Secondary particles propagate freely without any interaction during the forma- 
tion time r after the emission, 

E 

r = 7-ro = — -To, (A-17) 

m 

where the proper formation time tq is also one of the parameters of URASiMA. 
Values of the parameters a, (3 and tq are listed in Table V. 

A. 4. Comparisons with data of nucleus- nucleus collisions at high energies 

In order to examine the descriptive power of the event generator URASiMA-B, 
we made comparisons with the data of high-energy nuclear collision experiments. In 
this case, the initial state is given by the energetic nucleons forming the nuclei in free 
space. All parameters of the model are tuned to reproduce the data of hadron-hadron 
coUisions. ^^) " '^^^ 

The results of our simulations are compared with the experimental data of the 
E802 collaboration for the Si + Al central collision at 14.6 GeV/nucleon. ^^' Figure 
19 displays the rapidity distributions for the proton (left) and vr^ (right), where 
the squares denote the results of our simulation and the filled circles denote the 
experimental results. Figure 20 shows the rapidity dependence of the inverse slopes 
of transverse momentum distributions for the proton (left) and 7r+ (right). In both 
figures, URASiMA-B reproduces global features of experimental data quite well. 
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Fig. 19. A comparison between experimental results and our simulations. We plot the rapidity 
distributions of the proton (left) and tv (right) for Si + Al central collisions at 14.6 GeV/nucleon 
obtained in the E802 collaboration. Here, the squares denote the results of our simulation, and 
the filled circles denote the experimental results. 
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Fig. 20. A comparison between experimental results and our simulations. We plot inverse slopes 
of transverse momentum distributions of the proton (left) and n (right) for Si + Al central 
collisions at 14.6 GeV/nucIeon obtained in the E802 collaboration. Here, the squares denote 
the results of our simulation, and the filled circles denote the experimental results. 
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